﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace MathLib {
    public class DemogFilterAlg {

        private int periodCount;
        public bool fullind;
        private int groupCount;
        private GeneratedPopulation genPop;
        private double vnu, utf,utm;
        public Matrix P;
        public Matrix G;
        public Matrix AvgPopValue;
        public bool estind;

        public DemogFilterAlg(int periodCount, GeneratedPopulation genPop, double Vnu, double utf, double utm, bool fullind, bool estind) {
            this.periodCount = periodCount;
            this.genPop  = genPop;
            this.vnu = Vnu;
            this.utf = utf;
            this.utm = utm;
            this.fullind = fullind;
            this.estind = estind;
            this.groupCount = genPop.populations[0].NoRows;
            this.AvgPopValue = GenerationUtils.GenerateAvgPopValue(this.groupCount, this.genPop.populations);
        }




        public Matrix[] SolveWithFilter(Matrix[] Nu, Matrix[] MatrixNu) {
            Matrix[] x = new Matrix[periodCount];
            Matrix y;
            Matrix C = GenerationUtils.GenerateC(groupCount,fullind);
            Matrix Ut = GenerationUtils.GenerateUt(groupCount, utf, utm, AvgPopValue, genPop.populationMatrixes[0]);
            Matrix Vnu = GenerationUtils.GenerateVnu(groupCount, vnu, AvgPopValue, fullind);
            Matrix P = GenerationUtils.GenerateP(groupCount, 1);
            x[0] = genPop.populations[0];
            y = Matrix.Transpose(C) * genPop.populations[0] + Nu[0];
            Matrix A = genPop.populationMatrixes[0];

            for (int i = 1; i < periodCount; i++) {

                Matrix G = Matrix.Transpose(C) * P;
                G = G * C;
                G  = G + Vnu;
                if (Matrix.Det(G) != 0) 

                G = Matrix.Inverse(G);

                G = A * P * C * G;

                P = A * P * Matrix.Transpose(A) + Ut - G * Matrix.Transpose(C) * P * Matrix.Transpose(A);


                Matrix delta = y - Matrix.Transpose(C) * x[i - 1];
                Matrix newddelta;
                if (fullind)
                    newddelta = G * delta;
                else {
                    double s = delta[0, 0];
                    newddelta = G * s;
                }
               

                x[i] = A * x[i - 1] + newddelta;

                y = Matrix.Transpose(C) * genPop.populations[i] + Nu[i];

                A = genPop.populationMatrixes[i] + MatrixNu[i];
                if (estind) A = GenerationUtils.MortalityEstimation(x[i - 1], x[i],A, groupCount);

            }
            return x;
        }


        public Matrix[] Solve(Matrix[] Nu, Matrix[] MatrixNu) {
            Matrix y;

            Matrix[] x = new Matrix[periodCount];
            Matrix C = GenerationUtils.GenerateC(groupCount, true);

            x[0] = genPop.populations[0] + Nu[0];
            y = genPop.populations[0] + Nu[0];
            Matrix A = genPop.populationMatrixes[0] + MatrixNu[0];

            for (int i = 1; i < periodCount; i++) {

                x[i] = A * y;
                y = C * genPop.populations[i] + Nu[i];
                A = genPop.populationMatrixes[i] + MatrixNu[i];
            }
            return x;
        }
 
        public static Matrix CalclDeviation(Matrix m, Matrix[] estimate) {
            int minRow = estimate.Min(a => a.NoRows);
            int maxRow = estimate.Max(a => a.NoRows);
            if (minRow != maxRow)
                throw new Exception("Размерности векторов-оценок не совпадают");

            int rowCount = minRow;
            if (rowCount != m.NoRows)
                throw new Exception("Неверная размерность вектора мат. ожиданий");

            Matrix dev = new Matrix(rowCount, 1); ;
            for (int i = 0; i < rowCount; i++) {
                double d = 0;
                for (int j = 0; j < estimate.Length; j++) {

                    double s = m[i, 0] - estimate[j][i, 0];
                    d += s * s;
                }

                dev[i, 0] = Math.Sqrt(d / estimate.Length)/m[i,0];
            }
            return dev;
        }



        public static double CalcRelativelDeviation(Matrix m, Matrix estimate, Matrix disp) {

            double dev = 0;
            double sumReal = 0;
            double sumEst = 0;
          //  double sumNoise = 0;

            for (int i = 0; i < m.NoRows; i++) {
                sumReal += m[i, 0];
                sumEst += estimate[i, 0];
            }
            dev = Math.Abs(sumReal - sumEst) / sumReal;
            return dev;
        }


        public static Matrix CalcEffect(int expCount, int periodCount, GeneratedPopulation genPop, double Vnu, double Utf, double Utm, bool fullind, bool estind) {

            DemogFilterAlg df = new DemogFilterAlg(periodCount, genPop, Vnu, Utf, Utm, fullind,estind);

            Matrix m = new Matrix(2, 1);

            Matrix[] x = new Matrix[expCount];
            Matrix[] xx = new Matrix[expCount];

            double s1 = 0;
            double s2 = 0;

            int groupCount = x[0].NoRows;

            Matrix disp = new Matrix(groupCount, 1);

            disp = GenerationUtils.GenerateDisp(groupCount, df.AvgPopValue, Vnu);

            for (int i = 0; i < expCount; i++) {

                Matrix[] Nu = GenerationUtils.GenerateNoise(groupCount, genPop.populations, Vnu, fullind);
                Matrix[] MatrixNu = GenerationUtils.GenerateMatrixNoise(groupCount, genPop.populationMatrixes, Utf, Utm);

                x = df.Solve(Nu, MatrixNu);
                xx = df.SolveWithFilter(Nu, MatrixNu);

                double sigma1 = DemogFilterAlg.CalcRelativelDeviation(genPop.populations[periodCount - 1], x[periodCount - 1], disp);
                double sigma2 = DemogFilterAlg.CalcRelativelDeviation(genPop.populations[periodCount - 1], xx[periodCount - 1], disp);

                s1 += sigma1;
                s2 += sigma2;

            }
            s1 /= expCount;
            s2 /= expCount;
            m[0, 0]= s1;
            m[1, 0] = s2;
            return m;


        }
    }
}